###
### PARAMETERS 1 ###
###

# numero di veicoli
param K, integer, >= 0;

/* numeri di nodes (nodi meno deposito sdoppiato) */
param N, integer, >= 3;

# capacity
param C > 0;

# smallest ready time
param E >= 0;

#greatest due date
param L >= 0;

###
### SETS ###
###

# insieme dei veicoli
set vehicles := 1..K;

# nodes (V) e' l'insieme unione dei customers N unito al deposito {0, N+1}.
set nodes := 0..N+1;

# gli archi A sono sottoinsieme di VxV
set arcs within nodes cross nodes;

###
### PARAMETERS ###
###

# demand d[i] 
param d {nodes} >= 0;

# service time
param s {nodes} >= 0;

## time windows
# earliest service time
param a {nodes} >= 0;
# latest service time
param b {nodes} >= 0;

# coordinate dei customer per calcolo distanze (tempi)
param xcoord {nodes};
param ycoord {nodes};

# travel time t[i,j] fra vertice i e vertice j
# param t {arcs} >= 0;
param t {(i,j) in arcs} := sqrt((xcoord[i] - xcoord[j])^2 + (ycoord[i] - ycoord[j])^2); # Euclidean.

# bigM e' ha un valore elevato (ma il min possibile), tale da attivare  vincoli sulle finestre temporali qualora la variabile x_ij in soluzione sia unitaria 
param M {arcs} >= 0;

# non e' un parametro del modello matematico, ma serve solo per la visualizzazione del grafo con GraphvViz
param tripcolor {vehicles} symbolic;

param tsp_lb >= 0;
param aga_ub >= 0;

param graphFileName symbolic;

###
### VARIABLES ###
###

#  arc assignment
var x {arcs, vehicles} binary;

# arrival time
var w {nodes, vehicles} >= 0;

###
### OBJECTIVE ###
###

# eq 1.1, minimizzazione tempo totale
minimize totaltraveltime:
sum {k in vehicles, (i,j) in arcs} t[i,j] * x[i,j,k];

###
### CONSTRAINTS ###
###

# eq 1.2, da ogni nodo deve partire un veicolo
subject to assegnamento {i in nodes: i <> 0 and i <> N+1}:
sum {k in vehicles, (i,j) in arcs} x[i,j,k] = 1;

# eq 1.3, tutti i veicoli devono uscire dal deposito
subject to conservazioneflusso0 {k in vehicles}:
sum {(0,j) in arcs} x[0,j,k] = 1; #sum {(0,j) in arcs: j <> N+1 } x[0,j,k] = 1;


# eq 1.4
subject to conservazioneflusso{k in vehicles, j in nodes: j <> 0 and j <> N+1}:
sum{(i,j) in arcs} x[i,j,k] = sum{(j,i) in arcs} x[j,i,k];

# eq 1.5
subject to conservazioneflusson {k in vehicles}:
sum {(i,N+1) in arcs} x[i,N+1,k] = 1;

# eq 1.11 (sostituisce 1.6)
subject to linearizedtimeconstraints {k in vehicles, (i,j) in arcs}:
w[i,k]+s[i]+t[i,j]-w[j,k]<=(1-x[i,j,k])*100000;#M[i,j];

#eq 1.7a
subject to addresstimewindow1 {k in vehicles, i in nodes: i <> 0 and i <> N+1}:
a[i] * (sum {(i, j) in arcs } x[i,j,k]) <= w[i,k];

#eq 1.7b
subject to addresstimewindow2 {k in vehicles, i in nodes: i <> 0 and i <> N+1}:
w[i,k] <= b[i] * (sum {(i, j) in arcs } x[i,j,k]);

#eq 1.8a
subject to timeconstrainta {k in vehicles, i in nodes}:
E <= w[i,k];

#eq 1.8b
subject to timeconstraintb {k in vehicles, i in nodes}:
w[i,k] <= L;

#eq 1.9
subject to capacityconstraints {k in vehicles}:
sum {i in nodes} d[i] * (sum {(i,j) in arcs: i <> 0 and i <> N+1} x[i,j,k]) <= C;

# un upper bound da euristica
subject to ub:
	sum {k in vehicles, (i,j) in arcs} t[i,j] * x[i,j,k] <= aga_ub;

# un lower bound da
subject to lb:
	sum {k in vehicles, (i,j) in arcs} t[i,j] * x[i,j,k] >= tsp_lb;

solve;

## disegna grafo con GraphViz
printf 'digraph {\n' > graphFileName;
# i nodi
printf {i in nodes} '\t%d\t[label="%d\\n{%d} [%d,%d]" pos="%d,%d"];\n', i, i, s[i], a[i], b[i], xcoord[i], ycoord[i] >> graphFileName;
# gli archi
printf {k in vehicles, (i,j) in arcs: x[i,j,k] = 1 and (i <> 0 or j <> N+1)} "\t%d->%d\t[label=""%0.3f"" color=""%s""];\n", i, j, t[i, j], tripcolor[k] >> graphFileName;
#printf {(i,j) in arcs: (sum{k in vehicles} x[i,j,k]) = 0} "\t%d->%d\t[label=""%0.3f"" color=""gray""];\n", i, j, t[i, j] >> graphFileName;
printf '}\n' >> graphFileName;


printf {k in vehicles, (i,j) in arcs: x[i,j,k] = 1} "x[%d,%d,%d] = %f\n", i,j,k, x[i,j,k];
for {k in vehicles}
{
	printf "vehicle %d: 0 (%d)", k, w[0, k];
	printf {(i,j) in arcs: i <> 0 and x[i,j,k] = 1} ", %d (%f)", i, w[i,k];
	#param maxw := max {i in 1..N and x[i,j,k] = 1} w[i,k];    # maximum value of p
	#param minw := min {i in 1..N and x[i,j,k] = 1} w[i,k];    # minimum value of p
	#set S ordered by minw .. maxw;      # set up ordered set to hold p values
	#let S := setof {i in customer, } w[i, k];       # load set
	#let {i in I} p[i] := member(i, S);  # reload p 
	printf {(i,N+1) in arcs: x[i,N+1,k] = 1} ", %d (%f)\n", N+1, w[N+1,k];
}
printf "solution value: %f\n", totaltraveltime;

end;
